% Figure O4 (b)

axoptions={'scaled ticks = false',...
           'y tick label style={/pgf/number format/.cd, fixed, fixed zerofill,precision=2, set thousands separator={}}',...
           'x tick label style={/pgf/number format/.cd, fixed, fixed zerofill,precision=2, set thousands separator={}}',...
           'legend style={font=\normalsize}'}; %\scriptsize
load ../data/DataforMatlab/IP.csv
load ../data/DataforMatlab/torn.csv
load ../data/DataforMatlab\Uq_base.mat

Uq_base = Uq; 
%% Percentile

%% Percentile for the figures

WP_FS = -reshape(WP_FS_vec, [N, T]);
WP_total = -reshape(WP_total_vec, [N, T]);

Wfull = reshape(I_vec(:,:,:), [N, T]);
U = reshape(U_vec, [N, T]);
% Additional calculations
[~, indmax] = max(U);
[~, indmin] = min(U);
W = NaN(size(Wfull)); 
Wmax = Wfull(sub2ind(size(Wfull), indmax, 1:T));
Wmin = Wfull(sub2ind(size(Wfull), indmin, 1:T));
for tau = 1:T
    W(indmin(tau):indmax(tau), tau) = Wfull(indmin(tau):indmax(tau), tau);
end

logUmax = max(log(U));
logUmin = min(log(U));

PC = zeros(9, T);
for p = 1:9
    PC(p,:) = sum([Wfull <= IP(p,:)]);
end
PC(PC == 0) = 1;

Wq = zeros(9, T);
Uq = zeros(9, T);
WP_FSq = zeros(9, T);
WP_totalq = zeros(9, T);

for t = 1:T
    for p = 1:9
        Wq(p,t) = Wfull(PC(p,t), t);
        Uq(p,t) = U(PC(p,t), t);
        WP_FSq(p,t) = WP_FS(PC(p,t), t);
        WP_totalq(p,t) = WP_total(PC(p,t), t);
    end
end

Inf_base = log(Wq) - log(Uq_base);
Decomp1 = ((WP_totalq-WP_FSq))./Inf_base -1;
Decomp2 = ((WP_FSq))./Inf_base;
Decomp3 = (Inf_base - WP_totalq)./Inf_base;
Decomp_1 = [Decomp1(:,T)';Decomp2(:,T)']';
figure('name','Figure A4 (d)','NumberTitle','off')
bar(Decomp_1,'stacked')
legend('$\%\Delta$ inflation for X $\cup$ Y and X goods','Adj. due to $\Delta\log b_X$','Location','northwest','interpreter','latex' )
set(gca, 'XTick', 1:9, 'XTickLabels', {'10%','20%','30%','40%','50%','60%','70%','80%','90%'})
daspect auto
%matlab2tikz('../fig/Decomposition_const.tex','extraAxisOptions',axoptions);

    function F_X = interpNaN(X,Y)
        indexNaN = isnan(X);
        X(indexNaN) = [];% Eliminate NaN
        Y(indexNaN) = [];% Eliminate NaN
        [X, indexUN] = unique(X); 
        F_X = griddedInterpolant(X,Y(indexUN));
    end
